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We study the quantum dynamics of a number of model systems as their coupling constants are 
changed rapidly across a quantum critical point. The primary motivation is provided by the recent 
experiments of Greiner et al. (Nature 415, 39 (2002)) who studied the response of a Mott insulator 
of ultracold atoms in an optical lattice to a strong potential gradient. In a previous work, it had 
been argued that the resonant response observed at a critical potential gradient could be understood 
by proximity to an Ising quantum critical point describing the onset of density wave order. Here 
we obtain numerical results on the evolution of the density wave order as the potential gradient 
is scanned across the quantum critical point. This is supplemented by studies of the integrable 
quantum Ising spin chain in a transverse field, where we obtain exact results for the evolution of 
the Ising order correlations under a time-dependent transverse field. We also study the evolution of 
transverse superfiuid order in the three dimensional case. In all cases, the order parameter is best 
enhanced in the vicinity of the quantum critical point. 
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I. INTRODUCTION 

Recent experiments with ultracold atoms have 
achieved reversible tuning of bosonic atoms between 
superfiuid and Mott insulating states by varying the 
strength of periodic potential produced by standing laser 
light ^ . The physics of such ultracold atoms in the 
Mott insulating state can be described by bosonic Hub- 
bard model, well known in context of other condensed 
matter systems However, ultracold atoms in op- 

tical lattices offer much better control over microscopic 
parameters of the model. Consequently, it is possible 
to explore parameter regimes which are not available in 
other analogous condensed matter systems. 

This paper will focus on a particular experiment re- 
ported by Greiner et al. 0. With the boson system in 
the Mott insulating state, they applied a steep poten- 
tial gradient to the lattice, and observed its response. In 
a typical condensed matter system, one might have ex- 
pected a response analogous to that of a sliding charge 
density wave: no motion of atoms until a critical tilt was 
applied, and a sliding motion at all tilts above the critical 
tilt. However, the experiment observed strikingly differ- 
ent behavior: there was a strong resonant response in the 
vicinity of tilts where the potential energy drop between 
nearest neighbor optical lattice sites (E) equaled the re- 
pulsion between two atoms on the same site (U). For 
E ^ U, applying the tilt produced a noticeable change 
in the ground state, but (in contrast to sliding charge 
density wave systems), there was little change in the 
ground state for larger E until a second resonant peak at 
E ~ 2U. This resonant response is a clear indication that 
the atoms experience little extrinsic dissipation, and their 
dynamics should be described by an energy-conserving 
quantum Hamiltonian. 

A framework for describing the experiments of Greiner 
et al. was proposed in Ref . Isl (hereafter referred to as 
I). (We also note here the numerical studies of Braun- 
Munzinger et al. |^ which addressed these experiments 
by studying the time evolution of the underlying Bose- 



Hubbard model.) For w,\E — [/| ^ E,U, where w is 
the tunnelling matrix element between nearest neighbor 
lattice sites, it was argued that we need only focus on a 
set of states which were resonantly coupled to the origi- 
nal Mott insulating state. In one dimension, the resonant 
subspace could be described simply in terms of nearest 
neighbor dipole states, consisting of a particle and a hole 
excitation about the Mott insulator on nearest neighbor 
states; in higher-dimensions, the particle and hole were 
no longer constrained to be on nearest-neighbor sites but 
could reside anywhere on planes orthogonal to the po- 
tential gradient, but separated by a single lattice spac- 
ing. An effective Hamiltonian on such resonant subspaces 
was proposed in I, and its phase diagram was presented. 
In the regime of large potential gradient E — U > w, 
this effective Hamiltonian possessed ground states with 
density wave order with a period of 2 lattice spacings 
(see also Ref. for conditions under which other periods 
may obtain). It was argued in I that the proximity of 
the quantum critical point, associated with the onset of 
this density wave order, was responsible for the resonant 
response observed by Greiner et al.. 

The tilt experiments of Greiner et al. were carried out 
in highly non-equilibrium situations, and the approach 
of I was to describe these, to the extent possible, by an 
equilibrium analysis of an effective Hamiltonian describ- 
ing the primary states accessed over the experimental 
time scale. The purpose of the present paper is to di- 
rectly address the non-equilibrium dynamics of the tilted 
Mott insulator. We will mainly do this using the effective 
Hamiltonian of I. The specific question we shall address 
is the following. Begin with the system in the ground 
state in a regime of small E = Ei where there is no den- 
sity wave order. Then, suddenly change the value of E 
to a E — Ef, including values such that the ground state 
has density wave order at Ef. Allow the system to evolve 
under the resulting Hamiltonian. What is the nature of 
the state to which the system evolves at long times? We 
will find, as conjectured in I, that the density wave or- 
der that develops under this dynamic evolution is most 
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robust when Ef is near the quantum critical point. 

We will also address a similar question for the Ising 
chain in a transverse field, g. Like the models of L this 
model also has a regime, g < gc, where the ground state 
has spontaneous Ising order. However, this much simpler 
model is completely integrable, and so offers an oppor- 
tunity to analyze the non-equilibrium dynamics exactly. 
We initialize the Ising model in the ground state in a 
transverse field gi > gc- The transverse field is then 
changed rapidly to g = gf, and the wavefunction evolves 
at this field. We will compute equal-time correlations in 
this wavefunction as a function of the time t, including in 
the t ^ oo limit. In some cases, exact closed-form results 
will be obtained. The structure of these correlations as a 
function of gf bear some similarity to the results of the 
model of I as a function of Ef] however, there are some 
interesting differences which, we suspect, are related to 
the integrability of the Ising chain. 

We now outline the remainder of the paper. In Sec- 
tion^ we present numerical results on the dynamics of 
the one-dimensional dipole model of I. Section llTll will ad- 
dress the non-equilibrium dynamics of the Ising chain in 
a transverse field: this analysis uses the Jordan- Wigner 
transformation, and obtains the required dynamic cor- 
relation functions in the form of Toeplitz determinants. 
Section IIVI returns to the model of I, but turns to the 
dynamics in three dimensions; here, we use a combina- 
tion of mean-field theory and exact diagonalization to 
obtain results similar to those in Section m but with the 
order parameter now being a 'transverse superfluid' or- 
der. We review the results and discuss implications for 
experiments in Section Ivl 



II. DIPOLE DYNAMICS IN ONE DIMENSION 

This section will describe our numerical results on the 
quantum dynamics of the one-dimensional dipole model 
of the Mott insulator in a potential gradient. 

Starting from a parent Mott state with uq bosons per 
site, we identified the set of states which are resonantly 
coupled to the parent Mott state when U ^ E (recall 
that U is the repulsive energy between two bosons on the 
same site, and E, the 'electric field', is the potential drop 
between two nearest neighbor sites). In one dimension, 
the resonant subspace involves dipole states consisting 
of quasihole-quasiparticle pairs at adjacent sites, and the 
low energy behavior of the system can be described by 
the effective dipole Hamiltonian obtained in I: 

e 

+ {U - E)Y,d\d,. (1) 

The dipoles are subject to hardcore constraints that there 
is never more than a single dipole on any pair of nearest 



neighbor sites 

d\di < 1 ; dl^^dg+idlde — 0. (2) 

When the electric field E is adiabatically tuned through 
U, the ground state of the system changes from one with 
no dipoles {U S> E) to one with maximum possible num- 
ber of dipoles {E 3> U). At an intermediate critical 
electric field 

Ec = U+ IMOwy^noino + 1), (3) 

the system undergoes a quantum phase transition in the 
Ising universality class. 

As discussed in Section^ we study the dynamics of the 
ultracold atoms when the potential gradient is changed 
suddenly. Such a situation can be very easily achieved 
experimentally in these systems by rapidly shifting the 
center of the confining magnetic trap. We shall specifi- 
cally consider the situation where the change in the po- 
tential gradient is fast enough for the sudden perturba- 
tion assumption to be valid but slow enough to restrict 
the dynamics within the resonant subspaces so that the 
Hamiltonians 1^ (and H26() in Section irv|) are still valid. 

We assume that the atoms in the ID lattice are ini- 
tially in the ground state |\E'g) of the dipole Hamiltonian 
with E — Ei <^ Ec- This ground state corresponds 
to dipole vacuum. Consider shifting the center of the 
magnetic trap so that the new potential gradient is -E/. 
If this change is done suddenly, the system initially re- 
mains in the old ground state. The state of the system 
at time t is therefore given by 

exp{—ient / h)\n) , (4) 

n 

where \n) denotes the complete set of energy eigenstates 
of the Hamiltonian i?iD[-E/] in Q]), e„ = {n\Hi-D[Ef]\n) 
is the energy eigenvalue corresponding to state |n), and 
c„ = {n\'i'{t = 0)) = (nl^'c) denotes the overlap of the 
old ground state with the state \n). Notice that the state 
|\E'(t)) is no longer the ground state of the new Hamil- 
tonian. Furthermore, in the absence of any dissipative 
mechanism, which is the case for ultracold atoms in op- 
tical lattices, will never reach the ground state of 
the new Hamiltonian. Rather, in general, we expect the 
system to thermalize at long enough times, so that the 
correlations are similar to those of Hid [Ef] at some finite 
temperature. 

We are now in a position to study the dynamics of the 
Ising density wave order parameter 

e 

where N is the number of sites. The time evolution of O 
is given by 

0(i) = ^^C„^CnCOs[{E,n~ En)t/h] 
m.n 

x{m\Y,{-^Y4d,\n) (6) 
e 
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FIG. 1: Evolution of the Ising order parameter in © under 
the Hamiltonian //id[-E/] for no = 1. The initial state is 
the ground state of Hi-a[Ei\. All the plots in this section 
have U — 40, w — 1, and Ei = 32, and consequently the 
equilibrium quantum critical point is at Ec = 41.85. 



Eq.l^is solved numerically using exact diagonalization to 
obtain the eigenstates and eigenvalues of the Hamiltonian 
HiD[Ef]. Before resorting to numerics, it is however use- 
ful to discuss the behavior of 0{t) qualitatively. We note 
that ii Ef is close to Ei, the old ground state will have 
a large overlap with new one i.e. Cm ^ 5mi. Hence in 
this case we expect 0{t) to have small oscillations about 
0{t = 0). On the other hand, if Ef ;» E^ the two ground 
states will have very little overlap, and we again expect 
0{t) to have a small oscillation amplitude. This situation 
is in stark contrast with the adiabatic turning on of the 
potential gradient, where the systems always remain in 
the ground state of the new Hamiltonian HYD[Ef\, and 
therefore has a maximal value of (O) for Ej ^ Ec. In 
between, for Ej ^ E^, the ground state l^*) has a finite 
overlap with many states |m), and hence we expect 0{t) 
to display significant oscillations. Furthermore, if the 
symmetry between the two Ising ordered states is broken 
slightly (as is the case in our studies below), the time 
average value of 0{t) will be non-zero. 

This qualitative discussion is supported by numeri- 
cal calculations on finite size systems for system size 
N — 9, 11, 13. For numerical computations with finite 
systems, we choose systems with an odd number of sites 
and open boundary conditions, so that dipole formation 
on odd sites is favored, thus breaking the Z2 symmetry. 
The results are shown in Figs. El Fig. ^ shows the 
oscillation of the order parameter 0{t) for different val- 
ues of Ef for N = 13. In agreement with our qualitative 
expectations, the oscillations have maximum amplitude 
when Ef « 40 is near the critical value Ec — 41.85. For 
either Ef ^ Ec or Ef ^ Ec, the oscillations have a small 
amplitude around 0{t = 0). Furthermore, it is only for 
Ef f« Ec that the time-average value of 0{t) is apprecia- 
ble. Fig. 121 shows the system size dependence of the time 
evolution for Ef = U = 40. We find that the oscillations 



FIG. 2: System size (N) dependence of the results of Fig 
for Ef = 40. The curves are labelled by the value of A''. 
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FIG. 3: The curve labelled 'dynamic' is the long time limit, 
{0)t of the Ising order in ||SJ as a function of Ef (for N = 11), 
with other parameters as in Fig0 This long time limit can be 
obtained simply by setting m — n in JHJ. For comparison, in 
the curve labelled 'adiabatic', we show the expectation value 
of the Ising order O in the ground state of Hiu[Ef]; such 
an order would be observed if the value of E was changed 
adiabatically. Note that the dynamic curve has its maximal 
value near (but not exactly at) the equilibrium quantum crit- 
ical point Ec — 41.85, where the system is able to respond 
most easily to the change in value of E; this dynamic curve 
is our theory of the 'resonant' response in the experiments of 
Ref . Ill discussed in Section Q In contrast the adiabatic result 
increases monotonically with Ef into the E > Ec phase where 
the Ising symmetry is spontaneously broken. 



remain visible as we go to higher system sizes, although 
they do weaken somewhat. More significantly, the time 
average value of 0{t) remains non-zero, and has a weaker 
decrease with system size. In Fig.|3| we plot the long time 
limit of the Ising order parameter, (0)t, as a function of 
Ef, and compare it with the Oad, the value of the or- 
der parameter when E reaches Ef adiabatically and the 



4 



0.100 

<0>f 



0.075- 



0.050 



0.025 



0.000 




34 36 38 40 _ 42 44 46 



FIG. 4: Size dependence of the 'dynamic' results in Fig |3 
The sizes range from = 7 to A'^ = 15 (as labeled), with the 
intermediate values A*' = 9, 11, 13: {0)t decreases monotoni- 
caUy with A'' 



wavefunction is that of the ground state at E — Ef. We 
find that {0)t stays close to Oad as long as there is a large 
overlap with between the old and the new ground states. 
However, as we approach the adiabatic phase transition 
point, this overlap decreases and {0)t can not follow Oad 
any more. The deviation of {0)t is therefore a signature 
that the system is now in a different phase for the new 
value of the electric field. 

The 'dynamic' curve in Fig|31should be compared with 
Figs 5e,f in Ref. Q The latter show that the Mott in- 
sulator has a resonantly strong response to an applied 
potential gradient E ^ U . Here, we have found a sim- 
ilar resonant enhancement in a simple model system in 
one dimension, induced by the proximity of a quantum 
critical point. 

We comment briefly on the nature of the thermody- 
namic limit, iV — > oo for the results in Fig [301 For Oad 
it is clear that there is a non-zero limit only for E > Ec, 
when it equals the order parameter of the spontaneously 
broken Ising symmetry. If we assume that the system 
thermalizes at long times for the dynamic case, then (0)t 
corresponds to the expectation value of the equilibrium 
order parameter in Hio[Ef] at some finite temperature. 
In one dimension, it is not possible to break a discrete 
symmetry at finite temperatures, and so the thermody- 
namic limit of the order parameter must always vanish. 
By this reasoning, we expect (0)t to also vanish in the 
thermodynamic limit. This is consistent with the results 
in Fig. 21 where we show the N dependence of the long 
time limit (0)t. Our data are at present not extensive 
enough to definitely characterize the dependence of (O) < 
on N. 



III. DYNAMICS OF THE QUANTUM ISING 
CHAIN 



As a complement to the physically relevant, but numer- 
ical, computations in Section^ this section will describe 
similar results in a simpler, analytically tractable model. 
We will consider the integrable Ising chain in a transverse 
field, which also has a zero temperature, quantum phase 
transition between a phase with a broken Z2 symmetry 
and a symmetric phase. We will address questions on the 
evolution of the wavefunction under a time-dependent 
change in the transverse field. 

The model of interest in this section is 



(7) 



where uj'^ are Pauli matrices acting on a 'spin' on the 
sites, J, of an infinite chain. We have allowed the trans- 
verse field to acquire an arbitrary time dependence g{t). 
We will mainly consider here the case of a sudden change 
at time t — from an initial value .9(0^) = gt to a final 
value 5(0"'') — gf, but our methods easily generalize to 
the arbitrary time dependence in g{t). 

For time-independent g(t), Hj has a quantum critical 
point at g = gc — 1, with two equivalent ground states 
for g < gc related by a global Z2 spin- flip. However, un- 
like Section m we will not introduce any external pertur- 
bation which introduces a preference between these two 
states: all such perturbations destroy the integrability of 
Hi. Consequently, we do not obtain any useful informa- 
tion from the analog of the time-dependence of the order 
parameter in jSl, ©, as these quantities will be identi- 
cally zero at all times. Rather, we will compute here the 
two-point correlation function of the order parameter in 
an infinite chain, which is 

Gr.it) = {mw]'^]+n\m)- (8) 

Here \ip{t)) is the state of the system at time t, evolv- 
ing under the Schodinger equation specified by the time- 
dependent Hamiltonian Hj. In equilibrium, the infor- 
mation contained in a correlation function like is re- 
lated to an observable like that in © (which is the re- 
sponse in the Ising order parameter to perturbations in 
the boundary condition) by the fluctuation-dissipation 
theorem. However, we are not aware of any analog of 
such a theorem for the non-equilibrium case under con- 
sideration here, and so are not able to directly relate the 
results of the present section to those of Section ITU 

Our analysis of Hj proceeds with the standard Jordan- 
Wigner transformation, and we follow the notation and 
methods of Chapter 4 of Ref. 0. We express the S = 1/2 
states in terms of those of the spinless Jordan- Wigner 
fermion Cj, and after transforming to momentum space 
fermions Cfe, the Hamiltonian becomes 

Hi = ^ [2 (g - cos k) c\ck 

k 

-is\ak{c^ ^.c\ + c^kCkj " g ■ (9) 
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Now, transforming to the Heisenberg picture, we can fol- 
low the evolution of the system by solving the equations 
of motion 



^^f" • r 1 



(10) 



These equations are easily solved by a Bogoliubov trans- 
formation. Finally, the correlator in (|SJ) is computed by a 
simple generalization of the methods appropriate for the 
equilibrium case. A few details of such a computation 
appear in the Appendix. 

Here, we discuss the results for Gn (t) for the case of a 
sudden change from g{0~) = gi to 5(0+) — gj. For t < 0, 
we assume the system is in the ground state appropriate 
for g = gi^ and consequently G„(t < 0) is independent 
of t and equal to the well-known equilibrium result at 
g — gi- For t > Q, there is a non-trivial time dependence, 
and it is possible to obtain the general expression for 
Gn [t) as described in the Appendix. We will restrict our 
attention here to the simpler expression of the long time 
limit Gn{t ^00)7 which is the primary quantity of physi- 
cal interest. For this, we obtain the Toeplitz determinant 
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where 



with 



ar 
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2^ 



ikr ~ 



a(fc). 



(12) 



h{k) 



2(g/.9» + i)2-(g/+g»)(^' + i) 
- g/) 



[z ~ gi){\ - zg^Y 



(13) 



where z = e . 

We now need to evaluate the n x n Toeplitz determi- 
nant in Hll() . especially for the case of large n. In the 
equilibrium situation, this is aided by Szego's lemma, 
and its generalization in the Fisher-Hartwig formula 
For the present situation, the expression in l|13|) does not 
obey the winding number constraint required for applica- 
tion of the Fisher-Hartwig formula, and so we are unable 
to take advantage of this result. However, we shall show 
that an exact evaluation of (|11() is possible for two impor- 
tant special cases {gi = 00 and gi — 0), and supplement 
these by numerical evaluation of (|ll|l for other values of 

g«- 

In the case gi = 0, we have 



d{k) 



2z-gf{z^ + l) 
2{z-gf) 



(14) 



and it is straightforward to evaluate by contour inte- 
gration. This gives 
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(15) 



For the case gi — +co, we have 

2gfz - {z^ + 1) 



a{k) 



2(^-g/) 



(16) 



and Or is given by 





g/ < 1 


g/ > 1 


r < -1 
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^2 




r > 2 








(17) 



In both of these two cases, the following conditions are 
met: 

Condition 1 For gf>l,ar = for r < —1. 

Condition 2 For 5/ < 1, = for r > 2. 

Condition 3 For gf < 1, for r < —2. 

Using Condition 1, we can immediately write G'„(oo) = 
Qq for gf > 1. For gf < 1, define 



ao 



(18) 



so that Gn{oo) = 13°. Condition 2 gives DJ'j = 
a^rDj^_^ — aiDl^tX and Condition 3 gives = gfD^~^ 
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for r > 2. Also, = a_r- We can therefore write 



on 


_ ( ao 






\ 






_ 1 ao 






\a-i 


-ai9f 




_ I ao 


-fli 




\a-i 


-ai9f 



(19) 
(20) 
(21) 



This can be evaluated by diagonalizing the matrix. 

Collecting all the analytic results above, we have for 
the case gi = 0: 




n+l 


2" 



cosh 



G„(oo) = < 



1 



(n + 1) In 



1 



9f 



for 5/ < 1 

for gf > I 
(22) 

In the limit that n — > oo, the result for gf < 1 becomes 

ri+l 

G„(oo) ^ I I . (23) 



In the case gi = +oo, the corresponding results are 



G„(oo) 



cos (narccos(g/)) , for gf < I 



(25/) 



for 3/ > 1 



(24) 

Note that there are spatial oscillations in the correlator 
for the case where the field is reduced from a large posi- 
tive value [gi = +00) to a value below the critical point 
(5/ < l). 

Of these exact results, the case = 00 is the one 
that corresponds most closely to the physical situation 
discussed in Section ^] Here we start from an fully 'dis- 
ordered' initial state, and then suddenly change param- 
eters to values with increasing order (this is the analog 
of increasing E in Section^. For final parameter values 
gf > gc = ^, wc find here a result quite similar to that 
found in Section ^ from 124(1 we see that the order pa- 
rameter correlations decay with the a correlation length 
^f given by 







1 



ln(25/) ■ 



(25) 



This increases monotonically with decreasing gf, and is 
thus similar to the increase in the value of {0)t with in- 
creasing Ef for Ef < Ec in FigO By the analogy with 
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FIG. 5: Ising order correlations defined in ||HJ. The system 
is in the ground state of Hi for t < with g = gi — 2. At 
t — 0^ , the value of g is changed suddenly to g = gf, and 
remains at this value for all t > 0. Note that at long times, 
the order is best developed for gf = 1, which is the location 
of the equilibrium quantum critical point. This result is the 
analog of Figs |3] and 2] for the dipole model of Section ITTl 



Fig 121 we would expect here that there is a maximum in 
^f at g — gc- However, we find a somewhat different be- 
havior for gf < gc in 124|) : the correlations do not decay 
in a simple exponential, but now oscillate, with the pe- 
riod of oscillation becoming smaller with decreasing gf. 
So the correlations of the Ising ordered state are indeed 
best formed at = gc, but we find an unusual oscilla- 
tory decay of correlations for gf < gc. The oscillations 
are a clear indication of the absence of thermalization in 
the present model, and we expect they are special conse- 
quence of its integrability. 

We extended these analytic results by numerical eval- 
uation of (|ll|l for other values of gi, and found closely 
related behavior. Our results for gi = 2 are shown in 
Fig |S1 and these are the analog here of the results in 
Fig O and 0| As 5/ is decreased, the correlations be- 
come longer-ranged, until they reach a maximum range 
at gf = gc = 1. At smaller values ot gf, the correla- 
tions acquire an oscillatory behavior, but are also clearly 
shorter ranged. So the Ising order is best developed for 
gf near the quantum critical point. 



IV. DYNAMICS IN THREE DIMENSIONS 

We now return to the 'tilted' Mott insulator problem 
addressed in Section ^1 £^iid in I. Here we will address 
questions of quench dynamics for the three dimensional 
case. As discussed at length in I, the resonant subspace 
in 3D is described by quasiparticles and quasiholes which 
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are free to move in the directions transverse to the ap- 
pHed electric field. Consequently, the dipoles of Sec- 
tional which are bound quasihole-quasiparticle pairs in 
adjacent sites, constitute only a small part of the reso- 
nant subspace, and an effective Hamiltonian for unbound 
quasiparticle and quasihole states is necessary. A mean- 
field theory of this effective Hamiltonian was examined 
in I, and a fairly complex phase diagram was found. In 
addition to the Ising density wave order that appeared in 
one dimension, states with a transverse superfluid order 
were present. The latter states correspond to delocaliza- 
tion of the quasiholes and quasiparticles in the direction 
transverse to the applied electric field. 

In this section, we will address the quench dynamics 
across the transition associated with the onset of trans- 
verse superfluid order. This was found to be a second- 
order transition in the mean-field theory of I, and here 
we will extend the mean-field theory to an analysis of the 
non-equilibrium dynamics across the superfluid-insulator 
transition. We will not examine here the onset of Ising or- 
der, already studied in Section^ the present mean-field 
theory found a strong first-order transition for the on- 
set of Ising order. Our analysis will be restricted to the 
regime where both the superfluid and insulating states 
have no Ising density wave order. 

The effective mean-field Hamiltonian describing the 
dynamics of these quasiparticles and quasiholes can be 
written as in I: 



E 



wZ(nQhi(hi)* + (no -|- l)pe{pe) 



+h.cj - wy^no{nQ + l)(j>ehi-i + h.cj 

+ -^— — ^ (^plpi + h\ht^ - fii (^pl^-^^pi+l ~ h\he^ 

(26) 

Here ^ is a one-dimensional site index labelling sites along 
the longitudinal direction of the applied potential gradi- 
ent (the transverse degrees of freedom are treated in a 
mean-field approximation and so there is no dependence 
on the transverse site label), p and h are quasiparticle 
and quasihole annihilation operators, Z is the number of 
nearest neighbors in the transverse directions and /i£ de- 
notes chemical potential which enforces the constraints 



{p\+iPi+i) = {h\hi). 



(27) 



Although the Hamiltonian (|26|l has no non-linear terms, 
its diagonalization is non-trivial because of the hard-core 
constraint on all sites 



< 1, hW < 1, 



',h/ = 0. 



(28) 



The mean fields (pi) and (hi) correspond to transverse 
particle/hole superfluid order and were self consistently 



determined by diagonalizing the 3D Hamiltonian H26|) 
while maintaining l|28|l . 

We now consider the evolution of the ground state un- 
der a sudden shift in the value of E from E — Ei to 
E — Ef a.t t = . We place Ei in a regime where the 
ground state preserves all symmetry, and there is neither 
Ising or transverse superfluid order. The initial ground 
state I^E''^'^) will evolve according to the new Hamiltonian 
H3£,[{p),{h);Ef]. However, in contrast to the ID case, 
here the evolutions of the mean fields (p) and (h) have to 
be self-consistently determined. Within time-dependent 
Hartree approximation, we obtain 

,31)/ 



in 

dt 



y.{n\H:iYy[{pi{t)),{h,{t))-Ef\\m) 
(Pdt)) = Y.<nit)cnit){m\p\n) 

(heit)) = Y.''*rnit)cnit){m\h\n). 



(29) 



We used a basis of states |n) (the final results are, of 
course, independent of the choice of this basis) which 
are the complete set of eigenkets of the Hamiltonian 
HsuiiPi), {h{)',Ef], where (pj) and (/ij) are the ground 
states values of the particle and hole order condensates 
for E = Ef. All the states \n) maintain H28|) exactly, 
and so these hard-core constraints are fully respected by 
our calculation: this is what makes diagonalization of 
the Hamiltonian time consuming and numerically inten- 
sive. We note that these equations also maintain the 
constraints H27|l at all £ and t. 

We examined the above equations for the transverse 
superfiuid order using the same protocol used in Sec- 
tion ^1 for the Ising order. The set of Eqs. 1211 were solved 
self-consistently for longitudinal system size = 4. We 
consider the starting potential gradient Ei to be in the 
insulator phase with neither superfluid or Ising order, 
and ramp up the potential gradient to enter the super- 
fluid phase. The gauge symmetry of the superfluid order 
parameter is broken by adding a small symmetry break- 
ing term iJsym = — J2e v[{Pi + ^e) + where rj is an 
infinitesimal positive constant. In the absence of such 
a symmetry breaking term we would have (pe) and (hi) 
vanish identically for all t by gauge invariance. However, 
even an infinitesimal symmetry breaking is sufficient, in 
the suitable parameter regime, to induce appreciable val- 
ues of the superfluid order. In practice, this symmetry 
breaking is provided by the coupling of the system to its 
environment. 

With the symmetry breaking term present, the super- 
fluid order parameters (pe) and (hi) are initially real. As 
seen in Fig. El {h{t)) = ^f {hi{t)) develop coherent os- 
cillations once the potential gradient takes the value Ef. 
(These results are the analog of Figs I1I2I ') The oscilla- 
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FIG. 6: Oscillations of the hole superfluid order parameter. 
The plot is shown for f/ = 40, w = 1, no = 1, ^ = 4, = 20 
and Ef — 30. For these parameters, the quantum critical 
point is at Ec = 26.4. 
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FIG. 7: Analog of Fig but for the transverse superfluid 
order, using the same parameters (apart from Ef) as Fig 2] 
The "static" curve is the equilibrium superfluid order param- 
eter determined in I. The "dynamic" curve is the long time 
average of the real part of {h{t)). 



tions in {p{t)) are similar, but occur with a different pe- 
riod due to the inherent particle-hole asymmetry in Eq. 
1261 The long time average of the oscillations is purely 
real, and this is, of course determined by the symmetry 
breaking term. Such oscillations of the superfluid order 
parameter were also obtained receritly in a different con- 
text in Refs. Iglliol and by Levitov \u^ . 

We also examined the Ef dependence of the long time 
average of the superfluid order, and the analog of FigO 
appears in Fig \7\ Again the superfluid order is most 
strongly enhanced in the vicinity of the quantum criti- 



cal point. However, unlike Fig|31 we do not observe a 
precursor to the superfluid order in the insulating phase: 
this is surely an artifact of the mean-field treatment of 
the transverse degrees of freedom. 



V. CONCLUSIONS 

With advent of the study of quantum phase transitions 
in trapped atomic systems, there is a clear need for the- 
oretical studies in the highly non-equilibrium situations 
that experiments are often in. In particular, experiments 
can easily explore the change in the state of the system 
upon a sudden change in a parameter in the Hamiltonian. 
There are few general principles in such cases {e.g. there 
is no fluctuation-dissipation theorem which controls cor- 
relations of the final state), and theory is clearly still in 
its infancy. Two recent studies in this class |^ ^| , ex- 
amined the evolution of superfluid order under a sudden 
change in the optical lattice potential exerted on trapped 
bosons. 

It is clear that exact results on simple solvable models 
in non-equilibrium situations would be valuable. We have 
provided such an example here in Section IlIII where we 
examined the Ising chain in a transverse field, g. This 
model has a quantum critical point at 5 = gc, with 
spontaneous ferromagnetic order in the ground state for 
g < gc- We started the Ising model in the paramagnetic 
state {gi » gc), suddenly at t = changed g to a fi- 
nal value gf, and examined the long time development 
of correlations of the ferromagnetic order. (Our formal- 
ism also provided results for all t > 0, but we have not 
examined the detailed time evolution here.) The results 
are summarized in Fig [S| True long-range order does 
not develop at any value of gf, however, significant or- 
der parameter correlations do appear, and these are best 
formed for gj ~ gc- In a general non-integrable system we 
may expect thermalization at long times, at a tempera- 
ture such that the average energy equals that of the state 
at t = 0+. Such thermalization does not occur for the 
present integrable system, and the results have certain 
artifacts associated with this: the long-time correlations 
have an oscillatory spatial dependence for gf < gc. 

In the remainder of the paper we studied the non- 
equilibrium dynamics of models introduced in a previ- 
ous paper Q which addressed the response of a bosonic 
Mott insulator to a change in a strong potential gradient 
These models exhibit a number of quantum criti- 
cal points associated with the onset of Ising density wave 
and superfluid order. Our numerical studies here found a 
feature similar to that also obtained for the solvable Ising 
model: the order was best formed when the final param- 
eter value was in the vicinity of the associated quantum 
critical point, as illustrated in Fig|2| Here, and in Ref. 0, 
we have proposed this feature as the explanation for the 
resonant response observed by Greiner et al. Q upon 
'tilting' a Mott insulator of bosons in an optical lattice. 
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For t < 0, the field is gi and the system is taken to be 
in its ground state. We define the 7' fermions as those 
which diagonalize the Hamiltonian in the form l)A.l(l with 
this field. (Similarly, and FJ, are given by analogy with 
(TOll and (ESJ.) 

The state \tp) is therefore the vacuum of 7' fermions: 
in matrix notation, 



APPENDIX: COMPUTATIONS FOR THE ISING 
CHAIN 

The Jordan- Wigner transformation allows the Hamil- 
tonian of an Ising chain in a transverse field g to be writ- 
ten as 



(A.l) 



where 7fc is a fermionic annihilation operator (see Chap- 
ter 4 of Ref. 4). These are related to the Jordan- 
Wigner fermions Ck by a Bogoliubov transformation, 
parametrized by angle 9^, where 



tan 6k 



sin k 

g — cos k 



(A.2) 



In the present case, we define the 7 fermions as those 
that diagonalize the Hamiltonian for t > 0, with field gf. 

Since the Hamiltonian is throughout translationally 
symmetric, only fermionic states with opposite pseudo- 
momentum k and — k are mixed. We may therefore write 
the two component column vector 



(A.3) 



and similarly Ck for the Jordan- Wigner fermions. The 
Bogoliubov transformation relating Ck and F^. is ex- 
pressed as Ck = R^{9k)Tk, where 




R^(a) — cos 1- ia^ sin — 

^ ^ 2 2 



(A.4) 



and here tr^ is a 2 x 2 Pauli matrix. (These are used for 
conciseness of notation and should not be confused with 
the operators representing the 'spins' of the Ising chain.) 



■1 0' 
,0 0, 



2^ 



(A.5) 



Applying the Bogoliubov transformation to Ck and then 
Tk gives 

{^iTkTll^) = R'Hdk - 0'k)l{a^ + l)i?^(0, - 9'k). (A.6) 
Using K'"^ {a)a^ (a) = a^cosa — cr^ sina with 1/)^ ^ 
9k - 9'^ gives 



1 



(1 -I- cr^cos(/)fe - (7^sin(/)i;), (A. 7) 



as the set of matrix elements for the initial state. 

The time evolution of the operators now proceeds (us- 
ing the Heisenberg picture) according to the Hamilto- 
nian, (jXTJ, so that Ffc(t) = t/fc(i)Ffe(0), where 



Uk{t) 











R'H'^ekt). 



(A., 



The expectation values at any time can therefore be eval- 
uated using the algebra of SU{2) matrices. 

The n-site correlator can be written as (G„) = 
(Bo^i-Bi ■ • • -B„_iA„), where ^ c\ + q and = 
cl — Ci. Wick's theorem can then be used to write this ex- 
pression in terms of the expectation values of expressions 
bilinear in Ai and Bi. We therefore let 

= (^^^jjj^ = V2Ry (I) Ckit), (A.9) 

where Ak [Bk) is the Fourier transform oi Ai [Bi) so that 



UkA-k AkB^^k \ 1^^^ ^ {^\CkCl\^)W 




\BkA^k —BkB^k 
= 1 — cr^ (cos 9k COS (pk ~ sin (pk sin 9k cos 2ekt) 

; 9k cos 2ekt) + (sin (pk sin 2efei) 



hcr^ (sin 9k cos 0^ -f- sin 



sin (pk cos 9k I 



J sin 2efci 



, — e 



1 - sin(/)fe 
(cos (pk — i sin (pk cos 2ekt) 



—6'^*= (cos cj)k + i sin (pk cos 2ekt) 
1 + sin (pk sin 2ekt 



(A.IO) 

(A.ll) 
(A.12) 
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Transforming to real space and using the conservation 
of pseudomomentum gives 

{AiAj) = ^^e*'=('-^')(l-sin,^fcsin2efct) 
fc 

{BiBj) = i^^e^fe('-^)(-l-sin0fcSin2efet) 
k 

k 

+isin0fc cos2efet). (A. 13) 

The long time averages of these expressions are 



where 



M 
1 

2^ 



(A.15) 
(A.16) 



{BiB,) = -5ij 
{BiAj) = a;_j+i, 



(A.14) 



in the limit where the number of sites M becomes infinite. 
Here, a{k) = -e^(^^+'=) cos(6ifc -g;,). Wick's theorem then 
allows the expression for (G„) to be written as pi|) . 
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